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1. Introduction 

This paper presents some results on a well-known problem 
in Algebraic Signal Sampling and in other areas of applied 
mathematics: reconstruction of piecewise-smooth func- 
tions from their integral measurements (like moments, 
Fourier coefficients, Radon transform, etc.). Our results 
concern reconstruction (from the moments) of signals in 
two specific classes: linear combinations of shifts of a 
given function, and "piecewise Z?-finite functions" which 
satisfy on each continuity interval a linear differential 
equation with polynomial coefficients. 

Let us start with some general remarks and a conjecture. 
It is well known that the error in the best approximation of 
a C fc -function / by an iV-th degree Fourier polynomial is 
of order -^r. The same holds for algebraic polynomial 
approximation and for other basic approximation tools. 
However, for / with singularities, in particular, with dis- 
continuities, the error is much larger: its order is only . 
Considering the so-called Kolmogorow TV-width of fam- 
ilies of signals with moving discontinuities one can show 
that any linear approximation method provides the same 
order of error, if we do not fix a priori the discontinuities' 
position (see Q, Theorem 2.10). Another manifestation 
of the same problem is the "Gibbs effect" - a relatively 
strong oscillation of the approximating function near the 
discontinuities. Practically important signals usually do 
have discontinuities, so the above feature of linear rep- 
resentation methods presents a serious problem in signal 
reconstruction. In particular, it visibly appears near the 
edges of images compressed by JPEG, as well as in the 
noise and low resolution of the CT and MRI images. 
Recent non-linear reconstruction methods, in particular, 
Compressed Sensing ([ff] [3)) and Algebraic Sampling 
(||4] [12] [14] [6] [9]), address this problem in many cases. 
Both approaches appeal to an a priori information on the 
character of the signals to be reconstructed, assuming 
their "simplicity" in one or another sense. Compressed 
sensing assumes only a sparse representation in a cer- 
tain (wavelets) basis, and thus it presents a rather general 
and "universal" approach. Algebraic Sampling usually re- 
quires more specific a priori assumptions on the structure 
of the signals, but it promises a better reconstruction ac- 
curacy. In fact, we believe that ultimately the Algebraic 
Sampling approach has a potential to reconstruct "simple 
signals with singularities" as good as smooth ones. In par- 



ticular, the results of [TT] El El El strongly support 
(also apparently do not accurately formulate and prove) 
the following conjecture: 

There is a non-linear algebraic procedure reconstructing 
any signal in a class of piecewise C -functions (of one 
or several variables) from its first N Fourier coefficients, 
with the overall accuracy of order This includes the 
discontinuities' positions, as well as the smooth pieces 
over the continuity domains. 

At present there are many approaches available to a robust 
detection of discontinuities from Fourier data (see (8] 
[TP and references therein). The remaining problem seems 
to be an accurate estimate of the accuracy of the solution 
of the nonlinear systems arising. Our results below can be 
considered, in particular, as a step in this direction. On 
the other hand, they have been motivated by the results in 
IB El El, and in®|§|. 

2. Linear combinations of shifts of a given 
function 

Reconstruction of this class of signals from sampling has 
been described in J4]E). We study a rather similar prob- 
lem of reconstruction from the moments. Our method is 
based on the following approach: we construct convolu- 
tion kernels dual to the monomials. Applying these ker- 
nels, we get a Prony-type system of equations on the shifts 
and amplitudes. 

Let us restate a general reconstruction problem, as it ap- 
pears in our specific setting. We want to reconstruct sig- 
nals of the form 

N 
«'=1 3,1 

where the /j's are known functions of X = . . . , Xd), 
and the form ([TJ of the signal is known a priori. The pa- 
rameters dijj,, = (x° x , . . . , x J d ) are to be found from a 
finite number of "measurements", i.e. of linear (usually 
integral) functionals like polynomial moments, Fourier 
moments, shifted kernels, evaluation over some grid and 
more. 

In this paper we consider only linear combinations of 
shifts of one known function / (although the method 
of "convolution dual" can be extended to several shifted 
functions and their derivatives - see |[T6l ). First we con- 
sider general integral "measurements" and then restrict 



ourselves to the moments and Fourier coefficients. In what 
follows x = (xi, . . . ,Xd),t = (ti,. .. ,td), j is a scalar 
index, while k = (fci, . . . , kd), i = (ii, ■ ■ ■ , id) and n = 
(m, . . . , rid) are multi-indices. Partial ordering of multi- 
indices is given by k < k' k p < k' , p = 1, . . . , d. So 
we have 

s 

F(x)=J2 a jf( x + xj )' ( 2 ) 
3=1 

Let the measurements fik(F) be given by Hk{F) = 
J F(t)(pk(t)dt, for a certain (multi)-sequence of functions 
<p h (t), fc>0 = (0,...,0). 

Given / and ip = {(fk(t)}, k > we now try to find 
certain "triangular" linear combinations 



(3) 



0<i<fc 



forming, in a sense, some "/-convolution dual" functions 
(similar to a bi-orthogonal set of function) with respect to 
the system (fk(t). More accurately, we require that 



(4) 



We shall call a sequence ip = {V'fc(t)} satisfying (0, (0 
/ - convolution dual to ip. Below we find convolution dual 
systems to the usual and exponential monomials. 
We consider a general problem of finding convolution dual 
sequences to a given sequence of measurements as an im- 
portant step in the reconstruction problem. Notice that it 
can be generalized by dropping the requirement of a spe- 
cific representation (0: ipk(t) = Ej=o Ci,k<Pi{t). In- 
stead we can require only that J f(t)ipk(t) be expressible 
in terms of the measurements sequence fik- Also ipk in 
(0 can be replaced by another a priori chosen sequence 
r)k- This problem leads, in particular, to certain func- 
tional equations, satisfied by polynomials and exponents 
(as well as exponential polynomials and some kinds of el- 
liptic functions). 

Now we have the following result: 

Theorem 1. Let a sequence ij) = ipk (t) be f -convolution 
dual to ip. Define Mf. by Mfe = X)o<i<fc Q,kfM- Then the 
parameters a,j and x J in (0 satisfy the following system of 
equations ("generalized Prony system"): 



a-jipkix 3 ) = M fc , k = 0, .. 



(5) 



Proof We have M k = J2o<i<k Ct,kfJ-i = 
I Hi) Eo<i<fe C iik cpi(t)dt = jF(t)Mt) = 
Ej=i a i I f(t + x 3 )ip k (t)dt = ajifikix 3 )- 
In specific examples we can find the minimal number of 
equations in (0 necessary to uniquely reconstruct the pa- 
rameters a,j and in (0. 

2.1 Reconstruction from moments 

We are given a finite number of moments of a signal F as 
in (|2j in the form 



F(t)t n dt. 



(6) 



So here tp n (x) = x" 1 ■ ■ -x^ for each multi-index n = 
(ni , . . . , rid). We look for the dual functions ip n satisfying 
the convolution equation 



f(t + x)i> n (t)dt = x r 



(7) 



for each multi-index n. To solve this equation we ap- 
ply Fourier transform to both sides of ©. Assuming that 
f(w) G C°° (R d ), /(0) we find (see Q6)) that there is 
a unique solution to (O provided by 



lfin(x) = C, 



where 

Cn ,k 



1 



H7 







Qn— k 


1 


Q^n-k 


w=0 f(oj) 



(8) 



So we set the generalized polynomial moments as 

GnkfUk 



E 

k<7i 



(9) 



and obtain, as in Theorem[TJ the following system of equa- 
tions: 

s 

J^aji^)" ■ = M n , n > 0. (10) 

j'=i 

This system can be solved explicitly in a standard way 
(see, for example, lfl3l HI \l5l ). In one-dimensional case 
it goes as follows (see 1T31 ): from dTob we get that for 
z = (zi,...,Zd) the generalized moments generating 
function 



I(z) 



n£N d 



3=1 



d 1 

1=1 1 " 



xjzi 



(11) 



is a rational function. Hence its Taylor coefficients sat- 
isfy linear recurrence relation, which can be reconstructed 
through a linear system with the Hankel-type matrix 
formed by an appropriate number of the moments M n 's. 
This is, essentially, a procedure of the diagonal Pade ap- 
proximation for I(z) (see lfl3l ). The parameters a,j , x^ are 
finally reconstructed as the poles and the residues of I(z). 
For several variables the solution procedure is similar. 
In one dimensional case with the derivatives f"> included 
we have 



F(x) 



3=1 (=0 



(12) 



The corresponding moment-generating function in this 
case takes the form 



J w = EEE 

j=l 1=0 q=0 



l\ ( -l)^ au /{x^ 
(1 - xiz)i+ x 



(13) 



which is still a rational function (d-dimensional case with 
derivatives is similar). We would like to stress that in this 
case the dual polynomials ip k are not changed and they 
are given as in dU. Therefore also the formula for the 
generalized moments M n is the same as in (0. 



2.2 Fourier case 



In the same manner as in section 12.11 we now choose 
Indeed, 



e We get immediately tp k (x) = -jjTre lkx . 



f(t + x)Mt)dt= / f(t + x) 



1 



/(*) 



dt 



m 

/(*)' 



(14) 



Here the triangular system of equations © is actually not 

triangular any more but still since tpk{x) — -jy—<p-k(%) 

j (*) 

we can express the generalized moments through the orig- 
inal ones via Mk = -~j—p-k[F}. Now exactly as before 

/ (*■) 

we can find a generalized Prony system in the form 



1 



/(*) 



[F] = M k = J2 = J2 d5) 



The sequence {m*, = rn^g)} is given by the usual mo- 
ments 



m k {g) 



x k g{x)dx 



We subsequently formulate the following 

Piecewise D-flnite Reconstruction Problem. Given 
N,{ki\,K,,a,b and the moment sequence {rrik} of a 
piecewise D -finite function g, reconstruct all the parame- 
ters {ai t j},{£i},{aLi tn }. 

Below we state some results (see 0~| for detailed proofs) 
which provide explicit algebraic connections between the 
above parameters and the measurements {mk}. 
The first two theorems assume a single continuity interval 
(compare ITOll ). 

Theorem 2. Let K, = and D g = with 2) given by ( [TBI 
Then the moment sequence {mk(g)} satisfies a linear re- 
currence relation 



where pj = e~ lXj . In this case we get a rational expo- 
nential generating function and we can find its poles and 
residues on the unit complex circle as we did in the poly- 
nomial case. 



(E -a l) N (E -bl) N ■ E E ( fc < E ) ) m * = 

j=0 i=0 ' 

(18) 



2.3 Further extensions 

The approach above can be extended in the following di- 
rections: 1 . Reconstruction of signals built from several 
functions or with the addition of dilations also can be in- 
vestigated (a perturbation approach where the dilations are 
approximately 1 is studied in |fT31l ). 2. Further study of 
"convolution duality" can significantly extend the class of 
signals and measurements allowing for a closed - form sig- 
nal reconstruction. 

3. Reconstruction of piecewise Z> finite func- 
tions from moments 

Let g : [a,b] — > R consist of K. + 1 "pieces" go, .. .gic with 
K, > jump points 



a = £,o < f i . . . < & < 6c 



+i 



Furthermore, let g satisfy on each continuity interval some 
linear homogeneous differential equation with polynomial 
coefficients: £) g n = 0, n = 0, . . . , K, where 



N , kj 

s = E(E a ^)£j ( 

3=0 S=0 



dl 



(16) 



Each g n may be therefore written as a linear combina- 
tion of functions {ui]f =l which are a basis for the space 
AAs={/:S/ = 0}: 



A' 



g n (x) = ^2a hn u. l (x), n = 0, 1, (17) 

i=l 

We term such functions g "piecewise D-finite". Many 
real-world signals may be represented as piecewise In- 
finite functions, in particular: polynomials, trigonometric 
functions, algebraic functions. 



where E is the discrete forward shift operator and 
n^'^(fc,E) are monomials in E whose coefficients are 
polynomials in k: U^(k, E) = (-l)J ( /+t-j)! E '~ 3 - 

Theorem 3. Denote 



def 



E(E) = (E-aI) w (E-6I) JV , 



def 



= (f(E).n^(fc,E))m*, 



(0,j) k 



k=0 



Gj(x) = £(x)^—g(x) 



Assume the conditions of Theorem^ Then 

(1) The vector of the coefficients a = (ajj ) satisfies a lin- 
ear homogeneous system 



/ (0,0) 



Ha 



(1,0) 


(0,0) (1,0) 



(k N ,NY 
(k N ,N) 



\ 0,0 (1,0 
\ M 



M 



(fejv.A-) 
M 



( a 0,0 \ 



ai.o 



\ak N ,Nj 
(19) 



for all M G N. 

,(i,i) 



(2) v k l,J> = mj+fe (Gj(x)). Consequently, hj(z) is the 
moment generating function of Gj (x). 

(3) Denote Pj(x) =2^2iLo a i,j xl - Then the functions 
$ = {1, ho(z), . . . hjf(z)} are polynomially depen- 
dent: 

N 
3=0 

where Q(z) is a polynomial with degQ < max kj. 
The system of polynomials {z kj pj(z~ 1 )} is called the 
Pade-Hermite form for 



To handle the piecewise case, we represent the jump dis- 

, , def f X < 

continuities by the step function tl(x) = < and 

I 1 x > 

write g as a distribution 



g(x) =9o + ^2 9n(x)H(x - £ n ) 



(20) 



Theorem 4. Let K > and let g be as in (|20t with oper- 
ator D annihilating every piece y~^. Then the operator 



def 



n< 



(21) 



annihilates the entire g as a distribution. Consequently, 
conclusions of Theoretns\2\and\3\hold with D replaced by 
$) as in (ED- 

Proposition 5. Lef /C > a«if ^ e a basis for the 

space A/33, where D annihilates every piece of g. Assume 

x. 



dTTb a«(i denote c™ fc = " + 



1+1 x k u t (x) for n = 
A straightforward computation gives VM G N: 



\ 1,M 



-JV,0 



N.AI 



-N,0 



ajv,o 



\<xn,k) 



mi 
(22) 



The above results can be combined as follows to provide 
a solution of the Reconstruction Problem: 

(a) Let TV, {fcj}, /C, a, 6 and {nik(g)} be given. If /C > 0, 
replace S (still unknown) with !D according to (|2"TI) . 

(b) Build the matrix H as in ( fT9l ). Solve Ha = and 
obtain the operator D* = S) a which annihilates 5. 

(c) If JC > 0, factor out all the common roots of the poly- 
nomial coefficients of T>* with multiplicity N. These 
are the locations of the jump points {£ n }. The remain- 
ing part is the operator £> T which annihilates every g n . 

(d) By now and {£„} are known. So compute the basis 
for A/jjt and solve (l22l . 

The constants M and M determine the minimal required 
size of the corresponding linear systems ( TT9b and (l22l) in 
order for all the solutions of these systems to be also solu- 
tions of the original problem. It can be shown that: 

1 . There exists no uniform bound on M without any ad- 
ditional information on the nature of the solutions. 
Explicit bounds may be obtained for simple function 
classes such as piecewise polynomials of bounded 
degrees or real algebraic functions. 

2. For every specific D, an explicit bound M = M(D) 
may be computed for the system d22l . 

The above algorithm has been tested on exact reconstruc- 
tion of piecewise polynomials, piecewise sinusoids and ra- 
tional functions. 
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